namespace NBody.Analysis
{
    using System;
    using DataFile;
    
    public class DensityBin : Bin
    {
    	public override void Analyse(NBodySystem s)
        {
            switch (spacing)
            {
                case Spacing.Normal:
                    bin_normal(s);
                    break;
                case Spacing.Logarithmic:
                    bin_log(s);
                    break;
                case Spacing.EqualMass:
                    bin_equal_mass(s);
                    break;
            }
        }
        
        private void bin_log(NBodySystem s)
        {
            Console.Error.WriteLine("Binning density in radial shells with logarithmic spacing");
            
            double rmax = s.Rmax + 1e-10;
            double rmin = s.Rmin;
        
            double deltaR = Math.Log10(rmax / rmin) / (double)n;
        
            Console.Error.WriteLine("  Rmin = {0:f6}\n  Rmax = {1:f6}", rmin, rmax);
            
            data = new Table(n, 2, "Radius, Density");
            for (int i = 0; i < n; i++)
                data[i, 1] = 0.0;
            
            foreach (Particle p in s)
            {
                int pos = (int)((Math.Log10(p.Radius) - Math.Log10(rmin)) / deltaR);
                data[pos, 1] += p.Mass;
            }
            
            for (int i = 0; i < n; i++)
            {
                data[i, 0] = Math.Pow(10.0, Math.Log10(rmin) + 0.5 * deltaR + deltaR * i);
                double R = Math.Pow(10.0, deltaR * i) * rmin;
                double Rplus1 = Math.Pow(10.0, deltaR) * R;
                double volume = 4.0 / 3.0 * Math.PI * (Math.Pow(Rplus1, 3.0) - Math.Pow(R, 3.0));
                data[i, 1] /= volume;
            }
        }
        
        private void bin_normal(NBodySystem s)
        {
            Console.Error.WriteLine("Binning density in radial shells with normal spacing");
            
            double rmax = s.Rmax + 1e-10;
            double rmin = s.Rmin;
                                                                                
            double deltaR = (rmax - rmin) / (double)n;
            Console.Error.WriteLine("  Rmin = {0:f6}\n  Rmax = {1:f6}", rmin, rmax);                                                                                
            data = new Table(n, 2, "Radius, Density");
            data.Clear();
                                                                                
            foreach (Particle p in s)
            {
                int pos = (int)((p.Radius - rmin) / deltaR);
                data[pos, 1] += p.Mass;
            }
            
            for (int i = 0; i < n; i++)
            {
                data[i, 0] = rmin + 0.5 * deltaR + deltaR * i;
                double R = rmin + deltaR * i;
                double Rplus1 = rmin + deltaR * (i + 1);
                double volume = 4.0 / 3.0 * Math.PI * (Math.Pow(Rplus1, 3.0) - Math.Pow(R, 3.0));
                data[i, 1] /= volume;
            }
        }
        
        private void bin_equal_mass(NBodySystem s)
        {
            Console.Error.WriteLine("Binning density in radial shells with equal mass spacing");
            
            int num_per_bin = s.NumParts / n + 1;
            s.Sort();
            
            double rmax = s[s.NumParts - 1].Radius + 1e-10;
            double rmin = s[0].Radius;
            Console.Error.WriteLine("  Rmin = {0:f6}\n  Rmax = {1:f6}", rmin, rmax);
            
            data = new Table(n, 2, "Radius, Density");
            data.Clear();
            
            double m_enc = 0.0;
            double r_last = rmin;
            int c = 0;
            for (int i = 0; i < s.NumParts; i++)
            {
                m_enc += s[i].Mass;
                if ((i != 0 && i % num_per_bin == 0) || i == s.NumParts - 1)
                {
                    double r = s[i].Radius;
                    data[c, 0] = r_last + 0.5 * (r - r_last);
                    double volume = 4.0 / 3.0 * Math.PI * (Math.Pow(r, 3.0) - Math.Pow(r_last, 3.0));
                    data[c, 1] = m_enc / volume;
                    m_enc = 0.0;
                    r_last = r;
                    c++;
                }
            }
        }            
    }
}
